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lany  engineering  studies  require  the  evaluation  of  a Hilbert  Transform  of 
a tabular  function.  The  function  is  also  considered  as  an  even  function  in 
order  to  conform  to  physical  reality.  A numerical  procedure,  incorporated 
into  a Fortran  subroutine  HTRAN,  has  been  developed  to  accomplish  this  task. 
The  procedure  has  been  designed  to  avoid  the  theoretical  and  numerical 
singularities  and  employs  cubic  polynomial  interpolation  to  enhance  numerical 
accuracy. 
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The  Fortran  subroutine  HTRAN  computes  the  Hilbert  Transform  of  a tabular 
function  of  frequency  by  numerical  integration.  The  determination  of  the  Hilbert 
Transform  is  particularly  useful  in  engineering  applications  such  as  the  computa- 
tion of  the  complex  impedance  when  the  real  or  imaginary  part  is  known.  In  such 
instances,  the  real  and  imaginary  parts  are  related  by  the  Hilbert  Transform. 

In  the  physical  problems  under  consideration  we  can  consider  the  tabular 
function  of  frequency  as  arising  from  a Fourier  Transform  of  a causal  time  sys- 
tem, that  is,  a function  of  time  which  is  zero  for  t<  0.  Such  a system  gives  rise 
to  a complex  Fourier  Transform  with  a real  part  that  is  an  even  function  of  fre- 
quency and  an  imaginary  part  that  is  an  odd  function  of  frequency.  These  two 
parts  are  then  related  by  the  Hilbert  Transform.  Our  tabular  function  is  consid- 
ered the  even  function  of  frequency,  its  Hilbert  Transform  is  considered  the  odd 
function  of  frequency,  and  the  complex  combination  of  the  two  are  considered  the 
Fourier  Transform  of  a causal  system. 


! 
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2.  MATHEMATICAL  ANALYSIS 


For  a function  of  angular  frequency,  R(u),  the  Hilbert  Transform  is  defined 


X(w)  du' 

TT  J (0  - U 


The  integral  is  to  be  interpreted  as  the  Cauchy  Principal  Value  so  that  the  singu 
larity  arising  in  the  denominator  may  be  eliminated,  that  is 


X(w)  = 


4<sd_  + f -5<h!L 


du1  I . 


When  the  given  function  is  a function  of  frequency,  R(f),  we  have 


X(f)  = I J -j^rr  df 


which  is  interpreted  as 


X(f)  = 


f+e 


We  seek  to  evaluate  (3)  when  R(0  is  in  tabular  form  and  is  an  even  function  of 
frequency.  It  is  also  assumed  that  R(f)  is  described  tabularly  with  equal  spacing 
in  the  frequency  axis. 

Since  we  will  employ  numerical  Integration,  we  first  modify  expression  (3) 
in  order  to  eliminate  the  singularity.  We  note,  that 


f df 

J f - f ‘ 


i 


■■nuiriBi 


where  our  interest  in  f values  is  the  range  f^  < f £ Employing  the  evenness 
property,  we  are  able  to  obtain 


2f  r -R(n  ....  . r Rtf')  - Rtf)  ....  , r - Rtf) 

x(n  J fT77  df  J J UZJ d 


The  first  and  last  integrals  can  be  accomplished  directly.  Thus 


Xtf)  =-^-/n 


(1  -XV  fN 

V 1\V  1 t 2f  f _RJ 


f«)  - Rtf) 


df  . 


We  are  henceforth  interested  in  the  evaluation  of  the  integral  in  (12),  namely 


yhi  ■/  .SM^ajLd,. 


where  the  f values  are  in  the  range  f1  < f < f^-. 

3.  NUMERICAL  ANALYSIS 

Although  the  integrand  in  (13)  does  not  become  arbitrarily  large  anywhere, 
we  are  faced  with  an  indeterminacy  when  f equals  f.  The  integration  technique 
presently  described  has  the  feature  that  f1  never  attains  any  of  the  f values  at 
which  Y is  being  calculated,  thus  avoiding  the  indeterminacy  in  the  integrand. 

Consider  R(f)  defined  tabularly  as  pictured  by  the  solid  lines  in  Figure  1. 

(f.,  R.)  i = 1,  2, . . .,  N are  given  with  the  f values  equally  spaced.  The  subroutine 
computes  the  Hilbert  Transform  at  these  same  f values. 

Now  consider  another  set  of  f values,  denoted  by  f. , i * 1,  2, . ...  N-l,  each 
of  which  is  midway  between  two  f values.  Thus 


r fi  + fi+i 

i ~ 2 


= 1,  2 N-l 


as  pictured  by  the  broken  lines  in  Figure  1.  We  first  apply  (12)  to  obtain  the 
Hilbert  Transform,  X,  at  these  f values  where,  in  the  numerical  integration,  f 
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Figure  1.  Values  of  the  Independent  and  Dependent  Variables  Employed  in  the 
Numerical  Integration  Scheme 


ranges  over  the  f values.  Thus  the  indeterminancy  will  be  avoided  since  f will 
never  equal  one  of  the  f values. 

To  maintain  sufficient  accuracy  in  the  numerical  results,  we  employ  accurate 
cubic  interpolation  formulas  as  part  of  the  numerical  scheme. 

First,  we  need  to  obtain  the  values  R(f.)  i = 1,  2,  . . . , N-l  which  we  denote  as 
Rj  i = 1,  2,  . . . , N-l.  Employing  cubic  interpolation  for  the  interior  values  and 
parabolic  interpolation  for  the  end  values  we  obtain 

Rx  = (3R1  + 6 R2  - R3)/8 

Rt  = (-Rj.j  + 9R.  + 9Ri+1  - R.+2)/16  i = 2,  3 N-2 

Rjq_l  = (-Rjj_2  + ®^N-1  + * (15) 
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We  now  write  (13)  as 


fN 


r R(f')  - R(f;) 

Yi  = Y(fi)  = J " — ^ df' 


t,  r'2  - f? 


i = 1,2 N-l 


which  we  write  as,  employing  numerical  integration. 


N 


1 L I?  - 1 


i = 1,  2,  . . . , N-l 


j = l 3 


(16) 


(17) 


where,  when  N is  odd,  Simpson's  Rule  gives 
hl  r hN  = £f/3 


h.  = 4Af/3 

j = 2,4 N-l 

(j  even) 

(18) 

h.  = 2Af/3 

j = 3,  5,  . . . , N-2 

(j  odd) 

and  Af  is  the  spacing  in  the  frequency  axis.  1 

When  N is  even  the  Trapezoidal  Rule  is  used  to  include  the  last  interval. 

Thus 

h1  = Af/3 

h.  = 4Af/3 

j = 2,  4,...,  N-2 

(j  even) 

hj  = 2 Af/3 

j = 3,  5,  . . . , N-3 

(j  odd) 

(19) 

hN-1  = 5Af/G 


hN  = Af/2  . 

The  Hilbert  Transform  at  the  f.  i = 1,  2, ....  N-l  are  obtained  from  (12)  as 


xi  = x<V  =“T/n 


2f.Y. 

i i 


i = 1,2,...,  N-l 


(20) 
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To  obtain  the  Hilbert  Transform,  X.,  at  the  f.  i = 1,2,...,N  we  again  per- 
form an  accurate  interpolation;  cubic  interpolation  for  the  interior  values  and 
parabolic  interpolation  for  two  values  at  each  end.  Thus 


Xj  = (15XJ  - 10X2  + 3X3)/8 


X2  = (3X1  + 6X2  - X3)/8 


Xi  =<-Xi-2  + 9Xi-l  + 9Xi-Xi+l>/16 


i = 3,4, 


XN-la(“XN-3  + 6XN-2  + 3XN-l)/8 


XN  = (3XN-3  - 10XN-2  + 15XN-l)/8 


We  thus  obtain  the  tabular  results  (f.,  X^>  i = 1,  2, ....  N as  the  Hilbert 


Transform  values  for  the  given  (f.,  R^)  i = 1,  2, . . .,  N.  The  only  assumption 
made  concerning  the  function  R(f)  is  that  it  is  an  even  function  of  frequency, 
complying  with  phys;  al  reality.  It  is  also  assumed  that  the  function  R(f)  is 
described  tabularly  with  equal  spacing  in  f.  To  insure  sufficient  accuracy  in  the 
numerical  integrations  and  interpolations,  one  should  make  the  frequency  interval 
adequately  small  or,  correspondingly,  N sufficiently  large. 


4.  SUBROUTINE  DESCRIPTION 


The  user  employs  the  Hilbert  Transforms  subroutine  HTRAN  by  the  state- 


CALL  HTRAN  (R,  X,  N,  FBEG,  FEND) 


The  arguments  in  the  subroutine  are  described  as  follows: 

R - a dimensioned  array  containing  the  tabular  values  of  R. 

X - a dimensioned  array  containing  the  Hilbert  Transform  values  of  X 

returned  by  the  subroutine  HTRAN. 

N - the  number  of  values  contained  in  the  table  of  R (and  X). 

FBEG  - the  first  frequency,  fj,  for  the  R array. 

FEND  - the  last  frequency,  f^,  for  the  R array. 

Since  equal  spacing  in  the  frequency  axis  is  assumed,  the  items  N,  FBEG 


and  FEND  permit  the  subroutine  to  determine  the  values  fj  i = 1,  2, . . . , N. 


5.  PROGRAM  LISTING 


The  following  is  the  Fortran  listing  of  the  subroutine  HTRAN  as  written  for 
the  CDC  6600  at  Hanscom  Field,  Massachusetts. 

SUBROUTINE  HTRAN(R,  X,  N,  FBEG,  FEND) 

DIMENSION  R(3),  X(3) 

PI =3.  14159265359 
FDEL=(FEND-FBEG)/(N-1) 

F=FBEG+.  5*FDEL 
INC=MOD(N,  2) 

NI=N+INC-1 
NM1=N-1 
NIM2=NI-2 
DO  33  1 = 1,  NM1 
X(I)=0. 

IF  (I  .EQ.  1)  RX=(3.*R(l)+6.*R(2)-R(3))/8. 

IF  (I  .EQ.  NM1)  RX=(-R(N-2)+6.  *R(NM1)+3.*R  (N))/8. 

IF  (I  .EQ.  1 .OR.  I .EQ.  NM1)  GO  TO  20 
RX=(-R(I-l)+9.  *R(I)+9.  *Rd+D-Rd+2))/16. 

20  CONTINUE 
FI=FBEG 

DO  28  IP=1,  NIM2,  2 

X(I)=X(I)-r4.  *(R(IP+l)-RX)/((  FI+FDEL)**2-F**2) 

X +2.  *(R(IP  )-RX)/(  FI  **2-F**2) 

FI=FI+2.  *FDEL 
28  CONTINUE 
FEN=FEND 

IF(INC  . EQ.  0)  FEN =FEND-FDEL 
XU)=X(I)+(R(NI)-RX)/(FEN**2-F**2) 

X -<R<1  )-RX)/(FBEG**2-F**2) 

X(I)=FDEL/3.  *X(I) 

IF(INC  .EQ.  1)  GO  TO  30 

Xd)=Xd)+.  5*FDEL*((R(NI) -RX) / (FEN**2 -F!*t*2) 

X +(R(N)  -RX)/  (FEND®*2  -F*#2)) 

30  X(I)=2.  /PI#F*Xd)+RX/PI*ALOG 

X ((1.  -F/FEND)/(1.+F/FEND)*(F+FBEG)/(F-FBEG)) 

F=F+FDEL 
33  CONTINUE 
NM2=N-2 

XI=(15.  *X(1)-10.  *X(2)+3.  *X(3)))/8. 

X2=(3.  *X(l)+6.  *X(2)-X(3))/8. 

DO  31  1=3,  NM2 

XT=(-Xd-2)+9.*Xd-lD9.*Xd)-Xd+l))/16. 

Xd-2)=X1 

XI=X2 

X2=XT 

31  CONTINUE 

X(N)=(15,  *X(NM1)-10.  *XiNM2)+3.  *X(N-3))/8. 

X(N-1)=(3.  *X(NMl)+6.  *X(NM2)-X(N-3))/8. 

X(N-2)=X2 

X(N-3)=X1 

RETURN 

END 
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6.  NUMERICAL  EXAMPLES 

6.1  Example  1 

For  the  function 


R(f>  = 


1 0 < f < 1 

0 1 < f 


(22) 


with  R(f)  an  even  function,  the  exact  expression  for  the  Hilbert  Transform  in  the 
range  0 < f < 1 is  known  to  be 


X(f) 


1 - f 


1 + f 


0 <;  f < 1 


(23) 


This  example  reduces  to  the  trivial  case  since  the  integral  in  (12)  vanishes  and 
the  first  term  in  (12)  is  identical  with  (23). 

However,  exact  agreement  with  (23)  may  not  be  attained  by  the  subroutine 
since  interpolations  and  extrapolations  are  employed  to  find  the  Hilbert  Trans- 
form at  the  specified  f values.  Using  an  N value  of  501,  exact  agreement  was 
obtained  in  the  range  0 ' f < 1 except  for  values  of  f near  f = 1. 


6.2  Example  2 

For  the  function 


R(f)  = 


, l!2 

(i  - n o < r g i 

0 1 < f 


(24) 


with  R(f)  an  even  function,  the  exact  expression  for  the  Hilbert  Transfirm  in  the 
range  0 <C  f < 1 is  known  to  be 


X(f)  = -f  . 0 < f £ 1 


(25) 


Two  cases  were  examined  to  illustrate  the  accuracy  of  the  numerical  integra- 
tions. Table  1 depicts  the  comparison  of  the  two  cases  with  the  exact  results 
from  (25). 

We  first  note  an  overall  increase  in  accuracy  for  the  case  with  the  smaller 
spacing  along  the  f axis.  This  behavior  is  to  be  expected  with  a numerical  inte- 
gration scheme  such  as  the  one  employed. 
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Table  1.  Comparison  of  Example  2 With  Exact  Results  for  Two 
Different  Spacirgs  in  the  f Direction 


f 

exact 

Af  = .005 

Af  = . 002 

0 

0 

-.  9395268  X 10-10 

-.2351980  X 10“U 

0.  1 

-.  1 

-.  1000026 

-.  1000007 

0.2 

-.2 

-. 2000054 

-.2000014 

0.3 

-.3 

-.3000085 

-.3000022 

0.4 

-.4 

-.4000123 

-.4000031 

0.5 

5 

-.  5000172 

-.  5000044 

0.  6 

-.  6 

-.  6000242 

-. 6000061 

0.7 

-.7 

-.7000354 

-.7000090 

0.8 

-.8 

-.  8000572 

-.8000145 

0.9 

-.9 

-.9001214 

-.9000309 

1.0 

-1.0 

-1. 014396 

-1.009106 

We  also  note  a decrease  in  accuracy  in  both  cases  as  f increases  from  0 to  1, 
This  may  be  explained  by  the  nature  of  the  R(f)  function  (24).  The  magnitude  of 
the  slope  of  R(f),  and  consequently  also  for  the  integrand  in  (12),  increases 
greatly  as  f goes  from  0 to  1.  This  could  well  affect  the  accuracy  of  the  numer- 
ical integrations. 

6.3  Example  3 

For  the  function 


R(f)  = - 0 £ f (26) 

with  R(f)  an  even  function,  the  exact  expression  for  the  Hilbert  Transform  valid 
in  the  range  0 < f is  known  to  be 


cos  2 g f - 1 
' 2 * f 


0 < f 


(27) 


Two  cases  were  examined  here  to  illustrate  the  effect  of  the  numerical 
approximation  for  «.  In  one  case  the  upper  limit  of  f,  f^,  was  chosen  as  10  with 
an  N of  641  making  a spacing  Af  = 1/64.  In  the  second  case  fN  was  chosen  as  20 
with  an  N of  1281  thus  keeping  the  same  spacing  Af  = 1/64.  Table  2 shows  the 
comparison  of  the  two  cases  with  exact  results  (27)  for  the  upper  limit  of  f being  oo. 
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Table  2,  Comparison  of  Example  3 With  Exact  Results  for  Two  Different 
Approximations  for  f^  = » 


f 

exact 

o 

41 

J5 

fN=20 

0 

0 

. 7444819  X 10"4 

.7444819  X 10‘4 

0,  25 

-.  G366198 

6366234 

-.  6366199 

0.  50 

G366198 

-.  6366276 

-.  6366206 

0.  75 

2122066 

-.2122190 

-.  2122084 

1.0 

0 

-.  1647500  X 10'4 

-.2255452  X 10*5 

1.25 

-.  1273240 

-.  1273442 

-.  1273264 

1.  50 

21220G6 

-.2122311 

-.  2122094 

1.75 

09094568 

-.09097477 

-.  09094930 

2.00 

° 

-.3364389  X 10-4 

-.4225282  X 10-5 

As  was  to  be  expected,  the  case  with  the  larger  value  for  fM  produces  more 

i iN 

, accurate  results. 

i 

7.  SUMMARY 

The  subroutine  HTRAN  obtains  accurate  values  of  the  Hilbert  Transform  of 
a tabular  function  of  frequency,  which  is  equally  spaced  in  the  frequency  axis,  is 
an  even  function  of  frequency  and  is  zero  outside  its  range  of  tabular  definition. 

* Numerical  integrations  based  on  Simpsons  Rule  in  which  singularities  and 

j indeterminacies  are  eliminated  and  cubic  polynomial  interpolations  are  employed, 

j Values  of  the  Hilbert  Transform  are  obtained  for  the  same  frequency  values  as 

f are  specified  in  the  tabular  definition  of  the  function.  As  is  demonstrated  in  the 

j'  examples,  the  frequency  spacing  needs  to  be  small  in  order  that  the  numerical 

j integrations  and  interpolations  produce  accurate  results. 
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